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Abstract 

Several studies have shown that bursting neurons 
can encode information in the number of spikes per 
burst: As the stimulus varies, so does the length of 
individual bursts. The represented stimuli, how- 
ever, vary substantially among different sensory 
modalities and different neurons. The goal of this 
paper is to determine which kind of stimulus fea- 
tures can be encoded in burst length, and how those 
features depend on the mathematical properties of 
the underlying dynamical system. We show that 
the initiation and termination of each burst is trig- 
gered by specific stimulus features whose tempo- 
ral charactcristsics are determined by the types of 
bifurcations that initiate and terminate firing in 
each burst. As only a few bifurcations are possi- 
ble, only a restricted number of encoded features 
exists. Here we focus specifically on describing- 
parabolic, square-wave and elliptic bursters. We 
find that parabolic bursters, whose firing is initi- 
ated and terminated by saddle-node bifurcations, 
behave as prototypical integrators: Firing is trig- 
gered by depolarizing stimuli, and lasts for as long 
as excitation is prolonged. Elliptic bursters, con- 
trastingly, constitute prototypical resonators, since 
both the initiating and terminating bifurcations 
possess well-defined oscillation time scales. Fir- 
ing is therefore triggered by stimulus stretches of 
matching frequency and terminated by a phase- 
inversion in the oscillation. The behavior of square- 
wave bursters is somewhat intermediate, since they 
are triggered by a fold bifurcation of cycles of well- 



" Email: m.montcmurro@manchcstcr.ac.uk 



defined frequency but are terminated by a homo- 
clinic bifurcation lacking an oscillating time scale. 
These correspondences show that stimulus selec- 
tivity is determined by the type of bifurcations. 
By testing several neuron models, we also demon- 
strate that additional biological properties that do 
not modify the bifurcation structure play a minor 
role in stimulus encoding. Moreover, we show that 
burst-length variability (and thereby, the capacity 
to transmit information) depends on a trade-off be- 
tween the variance of the external signal driving the 
cell and the strength of the slow internal currents 
modulating bursts. Thus, our work explicitly links 
the computational properties of bursting neurons 
to the mathematical properties of the underlying 
dynamical systems. 



1 Introduction 

In the absence of noise, intrinsically burst- 
ing neurons generate periodic bursts. Peri- 
odic activity may serve to control basic phys- 
iological functions, as respiration (Smith et al., 
1991), digestion movements (Selverston and Miller, 
1980), egg laying (Alevizos et al., 1991), or in- 
sulin secretion (Meissner and Schmelz, 1974). It 
cannot, however, transmit information. As 
input noise increases, bursting responses be- 
come variable (Kuske and Baer, 2002; Su et al., 
2004; Pedersen and Srensen, 2006; Channell et al., 
2009), and the number of spikes per burst may fluc- 
tuate in time. Experimental studies performed in 
different sensory domains have demonstrated that 
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often, variability in burst length reflects variability 
in specific transient features of a time-dependent 
stimulus, for example, the orientation of a visual 
stimulus, (Martinez-Conde et ah, 2002), the loud- 
ness of an auditory stimulus, (Eyherabide et al., 
2009), or the velocity of a somatosensory stimulus 
(Arganda et al., 2007). Bursts of different lengths, 
therefore, are selective to different stimulus fea- 
tures. This paper aims at understanding which 
properties of bursting dynamics determine stimu- 
lus selectivity. 

The mapping between burst length and stim- 
ulus features is likely to depend on the bio- 
physical properties of the neuron, such as its 
morphology, the types of ionic channels lo- 
cated on the cellular membrane, their kinetics 
and their spatial distribution. These proper- 
ties may vary continuously in a complex, high- 
dimensional space. Yet, from the dynamical point 
of view, bursting neurons only display a finite 
number of qualitatively different behaviors, de- 
pending on the bifurcations governing the under- 
lying dynamical system (Izhikevich, 2007). Pre- 
vious studies (Izhikevich and Hoppensteadt, 2004; 
Mato and Samengo, 2008) have demonstrated that 
the computational properties of non-bursting cells, 
in particular their selectivity to specific stimulus 
features, are determined by the bifurcation govern- 
ing the onset of firing. We here hypothesize that a 
similar link between dynamical and computational 
properties can also be established for bursting neu- 
rons. 

By combining all the possible bifurcations that 
can initiate or terminate spiking, one can in 
principle construct 120 types of bursting neurons 
(Izhikevich, 2007). Most of these types, however, 
have never been observed experimentally, and some 
of them appear only rarely. We therefore focus on 
the three most ubiquitous types: parabolic, square- 
wave and elliptic. In order to explore the selectivity 
of each type, we simulate their activity when driven 
with a noisy input current and perform a statistical 
analysis of the stimulus stretches inducing bursting, 
obtaining the so-called preferred stimuli. To char- 
acterize the input-output mapping of each type, we 
search for the amplitude and frequency properties 
of the preferred stimuli and interpret those proper- 
ties in terms of the underlying bifurcations. 

In order to determine whether additional bio- 
logical details also shape stimulus selectivity, each 



type of burster is represented with three different 
models, capturing different amounts of biophysi- 
cal properties: a minimal normal-form model, a 
minimal conductance-based model, and a detailed 
conductance-based model. The main result is that 
stimulus selectivity is mainly shaped by the bifur- 
cation structure of the burster, and is roughly inde- 
pendent of further details. The present study thus 
links the dynamics and the computational charac- 
teristics of bursting neurons. 

Results 

Variability of bursting responses 

When driven with constant currents, intrinsically 
bursting neurons tend to generate periodic re- 
sponses. In the presence of noisy inputs, however, 
some variability in the number of spikes per burst 
may appear, as well as in the duration of the intra- 
burst and the inter-burst intervals. The amount of 
variability is not fixed, and example neurons can be 
found whose responses are remarkably periodic, or 
remarkably variable. We begin by identifying the 
factors determining the degree of periodicity 

The variables involved in the dynamical descrip- 
tion of bursting can be separated into two classes: 
fast variables and slow variables (Rinzel, 1987). 
The former participate in spike generation, and 
vary in time scales of the order of 1 ms. The latter 
vary at a much slower time scale (tens or hundreds 
of milliseconds). The evolution of the membrane 
potential is governed by equations of the form 

^ = /(^) + E / Uc + «E / bu rs t+^ x t, a) 

i 3 

where / is a function of the membrane potential 
(see Methods), /g piko and -^ urst represent the fast 
and slow currents respectively, I C xt is the external 
stimulus, and a is the coupling constant regulating 
the degree of coupling between fast and slow vari- 
ables. The evolution of the fast and slow currents is 
determined by additional variables, governed by ad- 
ditional equations whose functional form depends 
on the degree of biological detail (sec Methods). 

When fast and slow time-scales are well sepa- 
rated, slow variables define an almost constant cur- 
rent that competes with the external inputs synaps- 
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ing on the neuron. Throughout this paper, the ex- 
ternal current is modeled as 



4xt = h + cr£(t), 



(2) 



where 1$ is a constant DC, is a stochastic term 
of zero mean and unit variance (see Methods) and 
a represents the noise intensity. 

In Figure 1 we show that the amount of vari- 
ability in bursting responses is determined by a 
trade-off between the amount of noise a in the in- 
put signal and the strength a of the coupling be- 
tween fast and slow variables. One can therefore 
increase the response variability by cither increas- 
ing the amount of noise or decreasing the coupling 
between fast and slow variables. Panel (a) repre- 
sents the noiseless situation. The inter-spike inter- 
val (ISI) distribution contains two sharp peaks, the 
first one corresponding to the intra-burst intervals 
and the second one to the inter-burst silent peri- 
ods. The voltage trace looks perfectly regular. As 
we move to panels (b) and (c), the noise a increases. 
ISI distributions become wider, and voltage traces 
appear more disordered. If the amount of noise re- 
mains fixed, however, one can also vary the degree 
of disorder by modifying the coupling strength a. 
The larger the coupling, the larger the contribu- 
tion of the slow variables to the evolution of the 
membrane voltage. Noise therefore exerts a com- 
paratively weaker influence. Responses in panels 

(d) and (e) are calculated with the same a as in 
(b). However, increasing a (panel (d)) introduces 
order, resulting in a sharper ISI distribution and 
more regular spiking. Instead, decreasing a (panel 

(e) ) modifies the responses in a way that is simi- 
lar to the one obtained with increased noise (com- 
pare with panel (c)). Therefore, coupling strength 
and noise act as opposite forces in determining the 
amount of variability. 

Classification of bursting models in 
terms of their bifurcations 

In an adiabatic approximation, slow variables can 
be considered almost constant in short intervals. At 
the fast time scales, hence, slow variables operate as 
bifurcation parameters of the fast subsystem. One 
bifurcation initiates bursting and another one ter- 
minates it. The average values of the fast variables, 
in turn, self-consistcntly modulate the slow subsys- 
tem at longer time scales. The bifurcation structure 
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Figure 1: Periodicity vs. variability of burst- 
ing responses. A parabolic minimal conductance- 
based model (see Methods) is driven with a noisy 
input current of standard deviation a (expressed in 
fiA ms 1 ' 2 /cm 2 ) and mean value Iq. The magnitude 
of Jo was chosen to maintain the average firing rate 
at 75 Hz. (a) ISI distribution in logarithmic scale 
and example voltage trace for the noiseless case. 
The spike train is perfectly periodic, (b) and (c): 
increasing amounts of noise arc employed. The cou- 
pling strength is equal as in (a), (d) and (e): The 
amount of noise is fixed at 5 /iA ms 1 / 2 /cm 2 , but the 
coupling strength is increased in (d), and decreased 
in (e) , resulting in more ordered and disordered re- 
sponses, respectively. Noise and coupling strength 
thus act as opposite forces in determining response 
variability. 
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of different bursting neurons allowed Rinzel (1987) 
to classify bursters as parabolic, square- wave or el- 
liptic. 

In parabolic bursters, the initiation and termi- 
nation of spiking within each burst is governed by 
a saddle-node bifurcation on the invariant circle of 
the fast subsystem. These transitions are charac- 
terized by large-amplitude and low-frequency os- 
cillations. As a result, the first and last spike of 
the burst display ample voltage deflections. The 
ISIs, in turn, arc modulated throughout the du- 
ration of the burst, displaying a slow-fast-slow se- 
quence. This modulation gives parabolic bursts 
their name: The sequence of ISIs in a burst is a 
parabolic-shaped function of time. These features 
are displayed in Figure 2 (al), for the minimal 
conductance-based model. The minimal normal 
form and the detailed con-ductance-based models 
are topologically equivalent. Figure 2 (a2) shows 
the phase portrait of the parabolic burster pro- 
jected on the space defined by the two fast vari- 
ables, and one of the two slow variables. 

In square-wave bursters, the initiation of spik- 
ing is governed by a fold bifurcation, that is, a 
saddle-node bifurcation away from the invariant 
circle. This transition is characterized by large am- 
plitude and high-frequency spiking. Burst termina- 
tion occurs through a homoclinic bifurcation, dis- 
tinguished by large amplitude and low frequency 
fluctuations. Accordingly, the voltage trace in Fig- 
ure 2 (bl) shows a deceleration of spiking as the 
burst proceeds, without a significant attenuation 
of the height of voltage upstrokes. Figure 2 (b2) 
displays a burst in the 3-dimensional phase space. 

In elliptic bursters, the initiation of spiking is 
governed by two bifurcations, one of them respon- 
sible for the loss of stability of the fixed point rep- 
resenting the resting state, and the other for the 
creation of two limit cycles, the stable one repre- 
senting the spiking state, and an additional unsta- 
ble cycle. The fixed point looses stability through a 
subcritical Hopf bifurcation, and at approximately 
the same critical current value, the two limit cy- 
cles are created through a fold bifurcation of cycles. 
The combination of these two bifurcations is called 
a Bautin bifurcation (Arnold, 2008). Burst termi- 
nation involve the same two bifurcations occurring 
in the reverse direction: The fixed point becomes 
stable through a Hopf bifurcation and the limit cy- 
cles disappear through a fold bifurcation of cycles. 
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Figure 2: Properties of different bursting 
models. Minimal conductance-based models were 
stimulated with a constant external current, and 
the voltage trace is displayed (left) as well as a 
phase-space trajectory (right). Different types of 
bursters exhibit qualitatively different traces, (a): 
Parabolic burster, (al): Iq = 5^A/cm 2 . (a2): 
Io = 3/xA/cm 2 . (b): Square-wave burster, (bl) 
and (b2): I = 6/iA / cm 2 . C: Elliptic burster, 
(cl): I Q = 55.5/iA / cm 2 . (c2): I = b2^ik / cm 2 . 
Insets: amplification of the voltage trace in a single 
burst. 
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Elliptic bursts are characterized by large-amplitude 
and high-frequency fluctuations. As exemplified 
in the inset of Figure 2 (cl), the voltage trace 
exhibits subthreshold oscillations between bursts, 
characteristic of Hopf bifurcations. At spiking on- 
set, the voltage fluctuations suddenly increase in 
amplitude and switch to the frequency of the limit 
cycle created through the fold bifurcation of cycles. 
The subthreshold and suprathrcshold frequencies 
are associated with two different attractors, and 
therefore, do not necessarily coincide. The ter- 
mination process is qualitatively similar. Figure 2 
(c2) shows the typical subthreshold oscillations ob- 
served in elliptic bursters at the beginning and at 
the end of each burst appearing as 3-dimensional 
spirals. 

Bursting responses obtained with 
simple stimuli 

Before studying systematically bursting activity in 
response to stochastic inputs we briefly survey their 
behavior when driven with simple stimuli, as con- 
stant or sinusoidal input currents. 

When constant currents are employed (a = 0, in 
Equation 2) the number of spikes per burst typ- 
ically grows as a function of Iq, though not al- 
ways monotonically (Figure 3 (al, bl, cl)). In each 
case, Iq is varied between a minimum value, below 
which no spiking occurs, and a maximum value, 
beyond which spiking is no longer structured in 
bursts, i.e., tonic spiking begins. The mean intra- 
burst frequency remains roughly unchanged as Iq 
varies (panels (a2, b2, c2)). A constant mean intra- 
burst frequency, however, does not imply that all 
the ISIs inside the burst are equal. Parabolic and 
square-wave bursters, in fact, modulate their ISIs 
throughout the burst, so the standard deviation of 
the intra-burst frequency (shaded areas in Figure 3 
(a2, b2, c2)) is large, at least, when compared to 
elliptic bursters. The dented structure in panels 
(a2, b2, c2) corresponds to the steps in panels (al, 
bl, cl): When the number of spikes per burst in- 
creases in one unit, an extra spike is packed into 
each burst, so the intra-burst ISIs diminish, and 
so does the standard deviation. Finally, the inter- 
burst frequency grows with Iq (panels (a3, b3, c3)) 
with only small fluctuations each time the number 
of spikes per burst increases. 

When driven with subthreshold stimuli, 



parabolic and square-wave bursters do not 
display subthreshold oscillations, since their fixed 
points are foci. Elliptic bursters, instead, oscillate 
at the frequency associated with the elliptic fixed 
point (see Figure 2(c)). Therefore, when driven 
with subthreshold sinusoidal stimuli, elliptic 
bursters may resonate with the input frequency, 
as evidenced in the ZAP (impedance (Z) ampli- 
tude (A) profile (P), (Gimbarzevsky et al., 1984; 
Puil et al., 1986) depicted in Figure 3c. Parabolic 
(a4) and square-wave (b4) bursters have maximal 
voltage amplitude for zero-frequency stimulation, 
whereas elliptic bursters (c4) have a prominent 
resonance at 328 Hz. 



Classification of bursting models in 
terms of their degree of biological de- 
tail 

We now turn to the analysis of stimulus selectivity 
of different bursting models. Our working hypoth- 
esis is that the stimulus features associated with 
bursts of a given length depend on the dynamical 
properties of the bifurcations associated with the 
initiation and termination of bursts. Other bio- 
physical properties beyond the bifurcation should 
only weakly affect stimulus selectivity To test the 
hypothesis, we simulated each type of burster with 
three different models. The simplest models are for- 
mulated in terms of minimal normal-forms. They 
involve dimensionless variables, only loosely related 
to biophysical quantities. These models are de- 
signed to capture the topological properties of each 
bifurcation type in their purest form. As a first step 
towards biological realism we also used minimal 
conduct-ance-based models, formulated in terms 
of variables representing real biological quantities, 
though not intended to represent any specific neu- 
ron. These models were introduced by Izhikevich 
(2007) as the simplest conductance-based models 
with the bifurcations of real neurons. Finally, in 
an attempt to also include realistic neural models, 
we simulated detailed conductance-based models, 
proposed by Chay and Cook (1988), and further- 
analyzed by Bertram et al. (1995). These models 
represent pancreatic neurons that, depending on 
the value of the parameters, can exhibit parabolic, 
square-wave, or elliptic bursting. 
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Figure 3: Response characteristics of different types of bursting neurons. Minimal 
conductance-based models are stimulated with constant currents, and the number of spikes per burst 
(al, bl, cl) is measured as a function of the input strength. The intra (a2, b2, c2) and inter-burst (a3, 
b3, c3) frequencies are calculated as the inverse of the intra and inter-burst periods. The envelope of the 
voltage response to ZAP input currents is an estimation of the resonance curve for each model. Sub- 
threshold values of Iq have been used, near to the firing threshold. (a4): Iq = O.lfxA/ cm 2 , fluctuations 
of amplitude h = O.OlfxA/ cm 2 . (b4): I = 4/j.A/ cm 2 , h = 0.2^A/ cm 2 . (c4): I = 44fiA/ cm 2 , 
h = O.bfiA/ cm 2 . 



G 



Stimulus selectivity of different burst- 
ing models 

Stimulus selectivity was evaluated by calculating 
the burst-triggered average (BTA) of bursts con- 
taining a fixed number of spikes. To do so, we 
simulated several bursting models, differentiated 
by their type (parabolic, square-wave and elliptic) 
or by the amount of biophysical detail (minimal 
normal forms, minimal conductance-based mod- 
els, and detailed conductance-based models). In 
each case, bursts were identified by analyzing the 
ISIs in the spike train (see Methods). BTAs were 
calculated by pooling all the stimulus segments 
surrounding the generation of bursts of exactly n 
spikes and averaging them together. Therefore, 
each n-BTA represents the average stimulus around 
bursts of n spikes. We define t = as the time of 
burst initiation. 

In Figure 4, we show the results for the three 
parabolic models. Different curves correspond to 
bursts containing different number of spikes. The 
stimulus parameters are given in Tables 1-3. Al- 
though the detailed shape of the BTAs varies with 
the amount of biological detail, there are some gen- 
eral trends that are shared by all three models. In 
all cases, bursting is preceded by a marked increase 
in /cxt- In the minimal normal form and minimal 
conductance models (panels (a) and (b)), bursts 
are triggered by a peak that starts rising some 10 
ms before burst onset. The detailed conductance- 
based model (panel (c)) is governed by much slower 
time constants, and the peak begins to become no- 
ticeable 1-2 seconds before burst initiation, an in- 
terval 2 orders of magnitude longer than in the two 
preceding cases. Hence, the temporal scale of the 
structures in the BTA is determined by the details 
of the model, and not by the type of bifurcation. 

In the three models, for t < all curves overlap. 
No matter how many spikes are contained in each 
burst, bursts are always triggered by the same type 
of stimulus deflection. Hence, at the time of burst 
initiation, there is no way to predict how many 
spikes the burst will contain based on the shape 
of the stimulus thus far. To discriminate the stim- 
uli associated with short and long bursts, one must 
evaluate the shape of I ox t for positive times. Long 
bursts appear in response to sustained stimulation 
beyond burst initiation. Short bursts are termi- 
nated as soon as the stimulus drops. In the normal 
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Figure 4: BTA of parabolic bursters. Burst 
onset is at time zero, (a): Minimal normal form, 
(b): Minimal conductance-based model, (c): De- 
tailed conductance-based model. Error bars repre- 
sent the error of the mean. Different curves cor- 
respond to bursts containing different number of 
spikes. In the three models, bursts are triggered by 
a strong, depolarizing stimulus upstroke. The burst 
is prolonged as long as the excitation is maintained. 
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Figure 5: BTA of square- wave bursters. Burst Fi g ure 6: BTA of elliptic bursters. Burst on- 
onset is at time zero. Error bars represent the error set is at time zero. Error bars represent the er- 
of the mean. Different curves correspond to bursts ror of the mean - Different curves correspond to 
containing different number of spikes, (a): Mini- bursts containing different number of spikes, (a): 
mal normal form, (b): Minimal conductance-based Minimal normal form. Arrows show the position 
model, (c): Detailed conductance-based model. In of the upstroke that terminates the burst, (b): 
the three models, bursts are triggered by a small Minimal conductance-based model, (c): Detailed 
hyperpolarization, followed by a strong, depolariz- conductance-based model. In (b) and (c) arrows 
ing stimulus upstroke. To sustain the burst beyond show thc position of the downstrokc that termi- 
3 spikes, a second upstroke is required. nates thc burst - In thc thrcc models, bursts arc 

triggered by an oscillating input current of a well- 
defined frequency. Bursts are maintained until thc 
g input signal inverts its phase. 



form (panel (a)) and in the minimal conductance- 
based model (panel (b)), oscillations appear dur- 
ing sustained stimulation, and including the first 
large peak, there is one oscillation per spike in the 
burst. However, this characteristic is not a gen- 
eral property of all parabolic models, since it is not 
observed in the detailed conductance-based model 
(panel (c)). 

These characteristics may or may not appear 
in other types of bursting neurons. In Figure 5 
we show the BTAs corresponding to square-wave 
bursters. Also here bursts are triggered by a sharp, 
positive stimulus deflection, and once again, up 
to the time of burst onset (i < 0) the shape of 
the stimulus does not suffice to predict how long 
the burst will last. This time, we can also see a 
shallow hypcrpolarizing phase preceding the burst- 
triggering peak. As before, long bursts are asso- 
ciated with prolonged stimulation, whereas short 
bursts require a sharp downstroke to terminate. 
The shape of both the prolonged stimulation and 
the downstroke is remarkably similar across dif- 
ferent models (panels (a), (b), and (c)). How- 
ever, it differs from the shapes observed in the 
parabolic case (Figure 4). Long bursts in square- 
wave bursters require stimuli containing two sepa- 
rate peaks. In order to sustain spiking beyond the 
mean burst length (3 spikes per burst, for the mean 
currents used in Figure 5), an additional positive 
derivative in the current is required. This stimulus 
upstroke marks the initiation of the second peak 
in the BTA. In the minimal normal form, the sec- 
ond peak begins 0.5 ms after burst onset, and in 
the other two models, at 1 and 200 ms. In the 
three models, the timing between burst initiation 
and the beginning of the second peak coincides with 
the period of the spiking limit cycle at burst onset. 
Hence, to get long bursts we need a 2-peak oscilla- 
tion in the stimulus whose frequency coincides with 
the frequency of the limit cycle of the fold bifur- 
cation initiating spiking. The second peak in the 
stimulus therefore begins to be noticeable at the 
time when the second spike of the burst is fired. 

Elliptic bursting is markedly different from the 
previous two. In Figure 6 we see that the BTAs 
of the three elliptic models contain a pronounced 
oscillatory structure. The frequency of the oscil- 
lations remains roughly constant throughout the 
whole of the BTA in the normal form model ((a)). 
In the conductance-based models (panels (b) and 



(c)), the frequency diminishes after burst initiation. 
In all three cases, the frequency before burst initi- 
ation coincides with the frequency associated with 
the spiral fixed point, defined by the imaginary part 
of the eigenvalue loosing stability. The frequency 
after burst initiation coincides with the frequency 
of the spiking limit cycle. The two frequencies co- 
incide in the case of the normal form (a), but do 
not coincide in the conductance-based models (b) 
and (c). 

Burst initiation is triggered by gradually ampli- 
fying stimulus oscillations. Once again, the BTAs 
corresponding to bursts of different duration coin- 
cide for t < 0. Burst termination may be due to a 
marked upstroke (panel (a)) or downstroke (panels 
(b) and (c)), as indicated by the arrows in Fig- 
ure 6. In all cases, the terminating feature (be it 
an upstroke or a downstroke) is out of phase with 
respect to the ongoing oscillation. Hence, burst ter- 
mination is governed by a sudden phase inversion 
in the stimulus. The fact that the phase inversion 
may appear either as a peak or a trough implies 
that the extinction of the burst does not depend 
on the sign of the stimulus deflection, but rather, 
on the phase inversion. 

Dependence on simulation parame- 
ters 

Stimulus selectivity varies in a systematic fashion 
when the baseline current Jo, the amount of noise 
a, and the coupling strength a arc modified. In 
Figure 7 we display the BTAs of the normal-form 
models previously shown in Figures 4-6 but now for 
different coupling strengths a. 

In parabolic and square- wave bursters, as a 
changes, the family of BTAs obtained for different 
n values shifts either upward or downward. The 
direction of the shift is determined by the mean 
number of spikes per burst (n), which in turn, is a 
function of a. For a given a, the mean duration of 
each burst is (n). Bursts with n > (n) (panels (al) 
and (bl)) require sustained excitation in order to 
last longer than the typical length (n) . Bursts with 
n < (n) (panels (a3) and (b3)), instead, require a 
prompt downstroke to terminate firing rapidly, be- 
fore the length (n) is reached. In panels a2 and 
b2, (n) w 3, so some of the curves display down- 
strokes (those with n < (n)) and others sustained 
excitation (those with n > (n)). 
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Figure 7: BTAs for different coupling strengths a in normal-form models. Burst onset is 
at time zero. Error bars represent the error of the mean, sometimes too small to be visible. Different 
curves correspond to bursts containing different number of spikes ranging from n — 1 to n = 5 with the 
same color code as in Figures 4-6. (a) Parabolic bursting, al: a = 0.91. a2: a = 1. a3: a = 1.05. (b) 
Square-wave bursting, bl: a = 1.81. b2: a = 1. b3: a = 0.55. (c) Elliptic bursting, cl: a — 1.6. c2: 
a = 1. c3: a = 0.6. 
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The relationship between a and (n) is not uni- 
versal, i.e., different models show different trends. 
For example, increasing a results in a larger (n) 
in the square- wave model of Figure 7 (b), and a 
smaller (n) in the parabolic case of panel (a). Yet, 
in conductance-based models different trends may 
be found. The key factor determining whether (n) 
grows or diminishes with a is the sign of the av- 
erage value of the slow current. The current /burst 
exerts an influence on the fast sub-system that may 
be excitatory, inhibitory or neutral on average. In 
the parabolic model of Figure 7 (a) /burst is exci- 
tatory. In the square- wave model (panel b), /burst 
is inhibitory. Yet other models sharing the same 
bifurcations can be engineered with mean slow cur- 
rents of opposite sign. Therefore, the relationship 
between a and (n) is not a universal property of 
all parabolic or square-wave models. Once the re- 
lationship is given, however, stimulus selectivity 
varies systematically as a function of a. 

Elliptic bursters behave in the same way, though 
the effect is less visible by naked eye, because pro- 
longation and termination features are not easily 
separated visually. After burst onset, in-phase os- 
cillations constitute prolongation features. The last 
oscillation, appearing with inverted phase is a ter- 
minating feature. When n > (n), bursts need to 
be prolonged beyond their natural duration. The 
prolongation becomes increasingly difficult as time 
goes by, so the oscillations sustaining bursts during 
positive times increase in amplitude as the burst 
proceeds (marked with two lines in Figure 7 cl for 
n = 5). The terminating feature is instead small. 
When n < (n), bursts need to be terminated before 
they reach their natural duration. The termination 
signal must be particularly strong when the natural 
duration is large (large (n)). Therefore, the out-of- 
phase terminating feature in c3 (marked with ar- 
rows, one for each n) has a larger amplitude than 
in cl. 

Stimulus selectivity also varies when other pa- 
rameters are modified, as for example the stimulus 
mean Iq or the standard deviation a . In these cases, 
again, the critical factor is how these variations af- 
fect (n). Changes in (n) originated in changes in 
/o produce BTAs that, when arranged in increasing 
value of (n), are in all similar to the ones depicted 
in Figure 7. Changes in a also give rise to the 
same phenomena, with one additional effect: For 
large a, BTAs have less memory of events occur- 



ring in the distant past or distant future. As a 
consequence, BTAs look the same as in Figure 7, 
but with a damped envelope that gradually decays 
for positive and negative times. 

Independent components in stimulus 
selectivity 

The BTA is the average stimulus preceding a burst. 
When a single burst is generated, the triggering 
stimulus typically differs from the BTA up to a 
certain degree. Therefore, one can calculate not 
only the mean stimulus of each n, but also the 
variability around the mean. When studying a 
collection of scalar quantities, variability is cap- 
tured by the standard deviation. With vectorial 
quantities, a single number does not suffice to de- 
scribe variability, since there are many possible di- 
rections in which fluctuations may appear. Covari- 
ance analysis provides a tool for detecting the stim- 
ulus directions where uncorrclated variations are 
observed. Noticeably, variations are particularly 
large, or particularly small in only a few direc- 
tions. These are the so-called relevant directions: 
Stimulus directions where the standard deviation 
of the burst-triggering ensemble is either larger 
or smaller than in the prior stimulus (see Meth- 
ods). Covariance analysis also provides a system- 
atic procedure to reveal these relevant directions 
(Samengo and Gollisch, 2013). Here, the analysis 
is carried out for each intra-burst spike count n. 
The stimulus segments preceding bursts containing 
exactly n spikes are identified, and their temporal 
correlations are captured by the covariance matrix. 
The eigenvectors of this matrix represent the rel- 
evant features, and the corresponding eigenvalues 
provide a quantitative measure of the variance of 
the data in these directions (see Methods). Typi- 
cally, most burst-eliciting stimulus directions have 
approximately the same variance as the total stimu- 
lus, and only a few depart noticeably from the gen- 
eral trend. The eigenvectors associated with outlier 
eigenvalues constitute the relevant stimulus direc- 
tions. 

In Figure 8, we show the eigenvalues and eigen- 
vectors of the normal form of a parabolic burster. 
The two most prominent outlier eigenvalues are 
eigg and e2oo, and both of them have decreased 
variance. In panels (b-c) we see the eigenvectors 
associated with each of the outlier eigenvalues, for 
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bursts containing different numbers of spikes. 

Eigenvector V200 is the most significant relevant 
feature, since its eigenvalue has the largest devia- 
tion from unity. Its shape strongly resembles the 
BTA (Figure 4), except for the lack of a plateau at 
positive times, for long bursts. There is virtually no 
difference in the shape of the vectors correspond- 
ing to bursts containing different numbers of spikes. 
Variations in this direction capture the fluctuations 
in the stimulus deflection triggering bursts. Eigen- 
vector vi is therefore involved in shaping burst ini- 
tiation, but it contains no features related to burst 
termination. 

Parabolic bursting lasts for as long as the stimu- 
lus depolarizes the cell (Figure 4). That is, longer 
bursts are associated with more prolonged stimu- 
lation. Eigenvector V199 captures these differences: 
In Figure 8(b), we see that bursts containing a sin- 
gle spike are briefly stimulated, whereas 6-spikcs 
bursts are driven for almost 20 ms after burst on- 
set. Therefore, eigenvector Vigg crucially captures 
burst termination. 

Since burst initiation and burst termination are 
described by two different eigenvectors (V200 and 
Vigg, respectively), these two processes are uncorre- 
lated: Fluctuations in the features associated with 
burst initiation are not correlated with fluctuations 
responsible for burst termination. 

A similar situation is found for square- wave and 
elliptic bursters, as shown in Figures 9 and 10. In 
both cases, eigenvector V200 describes the shape of 
the stimulus features initiating the burst. In the 
case of a square-wave burster (Figure 9(c)), this 
feature is a depolarizing upstroke, whereas in the 
elliptic case (Figure 10(c)), bursting is initiated 
by an oscillation whose frequency coincides with 
the subthreshold resonance frequency of the cell. 
The remaining eigenvector (vigg in the square- wave 
burster, and Vi in the elliptic case) represents the 
stimulus features terminating the burst. 

When working with conductance-based models, 
similar results are obtained: Stimulus features re- 
sponsible for burst initiation and termination tend 
to appear in separate eigenvectors. The eigen- 
vectors responsible for burst initiation remain un- 
changed when n is varied. Those responsible for 
burst prolongation or termination, instead, mark 
the differences in stimulus selectivity of different n 
values. 
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Figure 8: Covariance analysis of a parabolic 
burster. Data obtained by simulation of the mini- 
mal normal- form model, (a): Example spectrum of 
eigenvalues for stimulus segments triggering bursts 
containing 5 spikes. The two most clear outliers 
(eigg and e2oo) are marked. Insets: eigenvector 
corresponding to each eigenvalue, (b-c): Eigenvec- 
tors associated with outlier eigenvalues, for bursts 
containing n spikes (one curve for each value of 
n). Burst onset is at time zero. Eigenvector V200 
only contains stimulus features initiating the burst, 
whereas vigg mainly describes structures needed to 
sustain and terminate the burst. 
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Figure 9: Covariance analysis of a square- 
wave burster. Data obtained by simulation of the 
minimal normal- form model, (a): Example spec- 
trum of eigenvalues for stimulus segments trigger- 
ing bursts containing 4 spikes. Two outliers (eigg 
and e2oo) are visible. Insets: eigenvector corre- 
sponding to each eigenvalue, (b-c): Eigenvectors 
associated with outlier eigenvalues for bursts con- 
taining n spikes (one curve for each value of n). 
Burst onset is at time zero. Eigenvector V200 con- 
tains stimulus features initiating the burst, whereas 
vigg describes the structures needed to sustain and 
terminate the burst. 
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Figure 10: Covariance analysis of an elliptic 
burster. Data obtained by simulation of the mini- 
mal normal-form model, (a): Example spectrum of 
eigenvalues, for stimulus segments triggering bursts 
containing 4 spikes. Two outliers (vi and V200) are 
visible. Insets: eigenvectors corresponding to each 
eigenvalue, (b, c): Eigenvectors associated with 
outlier eigenvalues, for bursts containing n spikes 
(one curve for each value of n). Burst onset is 
at time zero. Eigenvector V200 contains stimulus 
features initiating the burst, whereas vi describes 
the structures needed to sustain and terminate the 
burst. 



Discussion 



Initial studies on bursting neurons showed that 
they fire in a stereotyped fashion, and that 
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regular bursting appears spontaneously even in 
culture tissue (Chagnac Amitai and Connors, 
1989; Silva et al, 1991; Steriade, 1991; 
Agmon and Connors, 1992). Bursting was 
thus assumed to serve as a pacemaker signal in 
drowsy or anesthetized states (Steriade, 1991; 
Steriade et al., 1993), but could not encode tran- 
sient information, since a perfectly periodic signal 
has zero information rate. More recent studies, 
however, have shown that the number of spikes 
per burst may encode specific properties of the 
input stimulus. In sensory systems, these include 
the visual system of mammals (DeBusk et al., 
1997; Martinez-Conde et al., 2002), the auditory 
system of insects (Eyherabide et al., 2008, 2009; 
Marsat and Pollack, 2010), the somatosensory sys- 
tem of leech (Arganda et al., 2007) and the olfac- 
tory system of rodents (Cang and Isaacson, 2003). 
Theoretical studies have shown that the number of 
spikes per bursts encode the slope (Kepecs et al., 
2002) or the phase (Samengo and Montemurro, 
2010) of the input current driving a cell. 

Our analysis shows that these two behaviors con- 
stitute two extremes of a continuum phenomenon. 
Depending on the trade-off between the amount of 
noise in the external signal and the internal burst- 
related currents, neurons vary their coding capacity 
in a graded manner. When the external signal is 
constant, spike trains are perfectly regular. Neu- 
rons generate periodic bursts whose duration is de- 
termined by the oscillatory modulation produced 
by the slow variables. Since all bursts contain the 
same number of spikes, no information is encoded 
in burst duration. As fluctuations are incorporated 
to the external stimulus, bursts become irregular 
(see Figure 1). The degree of irregularity is deter- 
mined by a trade-off beetween the amount of noise 
in the input and the coupling between fast and slow 
variables. 

However, irregularity may or may not be in- 
formative. Bursts containing different number of 
spikes are only informative if they are selectively 
triggered by stimuli with specific properties. In 
order to assess whether there is a correspondence 
between burst length and stimulus attributes we 
calculated BTAs associated with different n values. 
We showed that the stimulus features associated 
with burst firing depend on the type of bifurcations 
initiating and terminating spiking in each burst, 
and are roughly independent of additional biologi- 



cal properties (Figures 4-6). 

The three explored models can be arranged along 
a resonator-integrator dimension. Elliptic bursters 
behave as emblematic resonators, as seen from 
their sharp subthreshold resonance and their pre- 
cisely tuned intra-burst frequency (Figure 2). Cor- 
respondingly, stimuli triggering, prolonging and 
terminating elliptic bursts have strong oscillatory 
components. Parabolic bursters behave diametri- 
cally opposite: They are archetypical integrators. 
Their two bifurcations have slow onsets and lack a 
characteristic time scale. Accordingly, stimuli trig- 
gering parabolic bursts are depolarizing currents 
that need to be sustained for as long as the burst 
lasts. Square- wave bursters lie somewhere in be- 
tween. The initiating bifurcation has a well defined 
frequency, but not the terminating one. Accord- 
ingly, stimuli triggering quadratic bursting have 
mixed characteristics: BTAs associated with long 
bursts contain two peaks, and the time interval be- 
tween burst onset and the initiation of the second 
peak coincides with the frequency of the firing limit 
cycle at burst onset. 

Covariance analysis showed that stimulus fea- 
tures triggering burst initiation appear in separate 
eigenvectors as those terminating bursts, implying 
that the two processes are uncorrelated. Moreover, 
stimulus features responsible for burst initiation do 
not depend on burst duration. At the time of burst 
onset, therefore, one cannot predict the duration of 
the burst based on the evolution of the stimulus 
thus far. Therefore, the generation of a burst - 
regardless of its length- reduces the uncertainty of 
the stimuli preceding burst generation, and thereby 
carries information about past stimuli. The length 
of the burst, in turn, reduces the uncertainty of 
the evolution of the stimulus during the duration 
of the burst, and thereby carries information about 
the events that take place between burst initiation 
and burst termination. 

Parameters as the DC stimulus component Iq, 
the standard deviation a and the coupling con- 
stant a affect the correspondence between n val- 
ues and stimulus attributes. The changes, how- 
ever, never mix different bursting models, that is, 
elliptic bursters are selective to oscillatory stimuli 
no matter the value of the parameters, and similar 
conclusions can be drawn for parabolic and square- 
wave bursters. When parameters are modified, the 
mean number of spikes per burst (n) varies. If 
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(n) increases, it becomes more difficult to obtain 
short bursts, and therefore, the few short bursts 
that appear occur in response to markedly pro- 
nounced terminating features. If several parame- 
ters are changed simultaneously in such a way as 
to leave (n) constant, stimulus selectivity remains 
roughly unchanged. 

Our results imply that in a given experiment, 
one may apply reverse correlation techniques to 
obtain BTAs and eigenvectors for different burst 
lengths, and from their shape, deduce the bifurca- 
tions governing the underlying dynamical system. 
Thereby, intrinsic properties of the recorded cells 
can be identified, even with extracellular measure- 
ments, where no information about the temporal 
evolution of the subthreshold activity is available. 
Care should be taken, however, to fulfill the follow- 
ing conditions. 

- The experiment should function in a range of 
mean resting potentials where the recorded cell 
bursts intrinsically. In some experiments, all 
observed bursts are driven by transient stim- 
ulus fluctuations, and the neuron stops burst- 
ing as soon as a constant stimulus is applied. 
Those cases are not described by the present 
theory, because the underlying dynamical sys- 
tem contains a single bifurcation: the one de- 
scribing the firing threshold. Such experiments 
should be analyzed with the theory developed 
for tonic neurons (Mato and Samengo, 2008). 

- The random stimulus should be appropriate 
for reverse correlation. In particular, it should 
have a high enough cutoff frequency, so that 
the missing high frequency components be ir- 
relevant to neuronal dynamics, that is, when 
present, they be filtered out by the capaci- 
tive properties of the cell membrane. When 
this condition is not fulfilled, BTAs contain os- 
cillations that do not represent intrinsic neu- 
ronal properties, but prominent stimulus fre- 
quencies, typically, the cutoff frequency. In 
order to avoid such artifacts, the required cor- 
rections for applying reverse-correlation tech- 
niques to colored stimuli should be made 
(Samengo and Gollisch, 2013). Even so, it 
should be noticed that those frequencies that 
are absent in the stimulus, are also absent in 
the obtained BTAs and eigenvectors, perhaps 



hindering the identification of prominent fre- 
quencies in the BTAs of elliptic bursters. 

- BTAs and eigenvectors should extend to posi- 
tive times long enough as to ensure that all ter- 
minating features are contained. They should 
also be discriminated by the number of spikes 
per burst so that terminating features can be 
seen gradually displaced in time, for longer 
bursts. 

- Noise should not be too large as to occlude the 
temporal evolution of BTAs and eigenvectors 
for positive and negative times. Recall that as 
noise increases, both BTAs and eigenvectors 
decay rapidly in time, so the features extend- 
ing to the past or the future may not be dis- 
cernible. Oscillations, for example, may damp 
out before a full cycle is completed. Since it 
may not be easy to decide, in a given exper- 
iment, which noise levels are appropriate, in 
doubtful situations, different runs using sev- 
eral noise levels are recommended. 

If these conditions are met, and assuming that 
the recorded cell produces either parabolic, square- 
wave or elliptic bursts, one may speculate which of 
the three alternatives best describes the data. El- 
liptic bursters are the best candidates for strongly 
oscillating BTAs with out-of-phase terminating fea- 
tures. Parabolic bursters constitute the best choice 
for purely depolarizing BTAs that last for as long 
as firing is sustained. Square-wave bursters are 
the appropriate choice when the upstroke triggering 
bursts is preceded by a shallow hyperpolarization 
and long bursts are sustained by a two-peak stim- 
ulus structure. 

Altogether, the function of bursters of the three 
major bifurcation types range from pure pacemak- 
ing to the encoding of information in the spike 
count per burst. In particular, burst onset and ter- 
mination are largely uncorrelatcd processes and can 
hence independently modify the trade-off between 
pacemaking and information transfer. Moreover, 
the initiating and terminating features are deter- 
mined by the specific bifurcations underlying burst 
initiation and termination. This implies that key 
computational aspects of bursting do not depend 
on the particular biological details of the neuron 
but rather on the basic topological structure of 
the underlying dynamical system. Physiologically, 
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bursters and the network they are embedded in 
have a wide repertoire of possibilities to regulate 
their function. In particular, our study suggests 
they can do so via changes in the coupling constant 
a, such as adaptation processes or short-term plas- 
ticity, as well as via changes to the effective input 
to the bursting neuron, such as the degree of syn- 
chronization in presynaptic neurons. 



Methods 



Burst identification 

Spike trains are parsed into sequences of bursts. 
Two consecutive spikes are assigned to the same 
burst or to different bursts, depending on their ISI. 
If the ISI is larger than a pre-defined threshold, the 
two spikes are assigned to different bursts. Oth- 
erwise, they are taken as part of the same burst. 
The threshold ISI is defined as the one where the 
ISI distribution reaches a local minimum (arrows, 
in the example of Figure 11 (a)). In all the simula- 
tions, the ISI distribution has a bimodal structure. 
This structure allows us to define the limiting ISI 
without ambiguity. Once spike trains are parsed 
into bursts, a distribution of n values is obtained 
(Figure 11 (b)). The mean and variance of this 
distribution depend on a, Iq and a. 

Stochastic stimulation 

In all cases, the external current is given by Equa- 
tion 2. There, is an Ornstein-Uhlenbeck pro- 
cess of zero mean and unit variance, with a correla- 
tion time r. In Figures 4-6 and 8-10, the value of Jo 
was chosen so that the average number of spikes per 
burst was approximately 3.5, and the value of a was 
chosen so that the standard deviation of the num- 
ber of spikes per burst was approximately 1.5. The 
value of t was chosen smaller than all the struc- 
ture in the BTAs, but not as small as to require 
too extensive simulations in order to obtain smooth 
BTAs. 
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Figure 11: Statistical properties of burst- 
ing neurons. Example simulation of the min- 
imal conductance-based models. Similar results 
are obtained with the normal forms, or the de- 
tailed conductance-based models, (a): ISI prob- 
ability densities, with the local minimum defining 
the threshold ISI marked by arrows. Smaller ISIs 
correspond to intra-burst intervals, and longer ISIs, 
to inter-burst intervals, (b): Probability of finding 
a burst containing n spikes, for different types of 
bursting neurons. 



Model bursting neurons 

For each type of bursting (parabolic, square-wave 
and elliptic), three neuron models are employed, 
with varying degree of biological detail. In all cases, 
with the parameters chosen for our simulations, a 
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positive stimulus input was required in order to ob- 
tain bursting responses. If the stimulus was suffi- 
ciently high, burst firing was replaced by tonic ac- 
tivity. 

Minimal normal-form models 

These models were designed to represent the rele- 
vant bifurcation in their canonical form. The equa- 
tions were taken from Izhikevich (2007), except for 
the dynamics of the slow current in the elliptic case, 
here modified to stabilize the subthreshold behav- 
ior, since in the original model, the slow current 
diverges when the fast currents are quiescent. 

The fast subsystem of parabolic and square- 
wave bursters is described by a single variable 
V evolving with the dynamics of a quadratic 
intcgrate-and-firc model neuron. This model con- 
stitutes the canonical form of any type I neuron 
(Ermentrout and Kopell, 1986). The equation is 



dV 
dt 



— = V Z + ah 



burst 



(3) 



where /burst is the current produced by the slow 
subsystem as detailed below, and I CKt is the exter- 
nal current. The coupling constant a is varied in 
Figure 7, and set equal to unity in all other plots. 
There is a resetting rule supplementing Eq. 3, stat- 
ing that whenever V increases above Vth, it is dis- 
continuously changed to a reset value Vr, and a 
spike is fired. 

The fast subsystem of elliptic bursters involves 
two variables. The equations are usually given in 
radial coordinates (r, #) , such that the variable rep- 
resenting the voltage is rcosi?. The dynamics is 
governed by 



dr 

dt 
dfi 

dt 



/burst r + cr 3 + dr 5 
1 



(4) 



The system of Eqs. 4 is the canonical form of a type 
II neuron with a subcritical bifurcation at firing on- 
set (Brown et al., 2004). The external current Text 
is incorporated as an additive term in the equation 
for the voltage, obtained by transforming Eq. 4 to 
cartesian coordinates. 

We now describe the dynamics of the slow sub- 
systems. For parabolic bursters, /burst = Ui — U2, 



where 

d«i r/ , 

— — = -/illil + did{t — Espikejj 

^ = -(i2U2 + d 2 5(t-t spikc ), (5) 

where i S pike represents the time when the V crosses 
the threshold voltage Vth and generates a spike. 
At those times, the slow variables u\ and u 2 are 
discontinuously reset to u\ + d\ and u 2 + d 2 , re- 
spectively. For square- wave bursters, /burst = — u\, 
where 



du\ 

— = -/llli! + dxO\t - t S pikc) 



(6) 



When V crosses the threshold Vth, the slow vari- 
able u\ is discontinuously changed to u\ + d\. For 
elliptic bursters, 



d/bu 



dt 



"Ml (/burst + Ar 2 ) 



(7) 



In the case of elliptic bursting, no discontinuous re- 
setting rule is used. Instead, a spike is fired when- 
ever r cos $ crosses the threshold 0.75. Table 1 spec- 
ifics the parameters used in the simulations. 



Parameter 


Parabolic 


Square- wave 


Elliptic 




bursting 


bursting 


bursting 


Io 


-0.1 


-0.1 





a 


0.25 ms 1 ' 2 


1.4 ms 1 / 2 


0.75 ms 1 / 2 


T 


1 ms 


0.5 ms 


0.2 ms 


Vth 


20 


10 




Vr 


-1 


1 




Mi 


0.1 ms -1 


0.015 ms- 1 


0.0025 ms" 1 


di 


1.1 


0.22 






0.02 ms" 1 






d 2 


0.55 






A 






1.25 


c 






0.4 


d 






-0.2 


dt 


0.1 ms 


0.05 ms 


0.01 ms 


a 


1 


1 


1 



Table 1: Parameters of the minimal normal-form 
models used in Figs. 4-6 and 8-10. Figure 7 was 
built with the same parameters except for the cou- 
pling strength a. 
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Minimal conductance-based models 

These models are taken from Izhikevich (2007). 
Spiking is induced by the fast subsystem, contain- 
ing a fast persistent sodium current and a potas- 
sium current. Bursting is modulated by the slow 
subsystem, involving one or two additional slow 
variables. For all bursting types, the equations of 
the fast subsystem read 



r dV 



An 



"ffNa m^V) {V - VNa) - ffK n (V 
-ffh (V - Vl) + Ct/burst + 4xt 

ioo{V) - n] /r n . 



V K ] 



In Figure 1, a takes the values 0.5, 1, and 1.1. 
Throughout the rest of the paper, a is maintained 
fixed at 1. The activation curves Xoo are 



i 



exp 



(V? /2 - V)/k* 



(9) 



and the parameters depend on the type of bursting, 
as specified in Table 2. In the parabolic burster, 
the slow subsystem is described by two variables, 
associated with sodium and potassium currents, re- 
spectively That is, the bursting current reads 

4urst = -gm S m s (V - VNa) - 0ns ^(V - Vk), (10) 



and 



dm s 
dn s 



[m s ooiV) - m s ] /r m 
[n s oo(V) - raj /r ns . 



(11) 



In the square-wave and the elliptic bursters, the 
slow subsystem is described by a single variable, 
associated with a potassium current, 



^burst = -9n S n s (V - Vk), 



(12) 



and n s is governed by Eq. 11. Parameters are given 
in Table 2. 

Detailed conductance-based models 

The Chay-Cook model is a detailed conductance- 
based model of a pancreatic /3-cell that can produce 
parabolic, square- wave or elliptic bursting, depend- 
ing on the parameters (Chay and Cook, 1988). The 



Pcirfiinctcr 


P&r&bolic 


ScUlclXG-WclVG 


Elliptic 




bursting 


bursting 


bursting 


I 


1.2 fiA/cm 2 


4 /iA/cm 2 


48 ^A/cm 2 


a 


2 ^A/cm 2 


2.75 iik/cm 2 


3 ilk/ cm 2 


T 


1 ms 


0.5 ms 


1 ms 


c 


1 /iF/crn 2 


1 /iF/cm 2 


1 /iF/cm 2 


VNa 


60 mV 


60 mV 


60 mV 


Vk 


-90 mV 


-90 mV 


-90 mV 


v L 


-80 mV 


-80 mV 


-80 mV 


5Na 


20 mS/cm 2 


20 mS/cm 2 


4 mS/cm 2 


9K 


10 mS/cm 2 


9 mS/cm 2 


4mS/cm 2 


9h 


8 mS/cm 2 


8 mS/cm 2 


1 mS/cm 2 


V l/2 


-20 mV 


-20 mV 


-30 mV 


k m 


15 mV 


15 mV 


7 mV 


v n 

V l/2 


-25 mV 


-25 mV 


-45 mV 


k u 


5 mV 


5 mV 


5 mV 


T n 


1 ms 


0.152 ms 


1 ms 


9ms 


3 mS/cm 






yms 
V l/2 


-40 mV 






f,ms 


5 mV 






'ms 


20 ms 






5ns 


20 mS/cm 2 


5 mS/cm 2 


1.5 mS/cm 2 


V l/2 


-20 mV 


-20 mV 


-20 mV 


k ns 


5 mV 


5 mV 


5 mV 


7"ns 


50 ms 


20 ms 


60 ms 



Table 2: Parameters of the minimal conductance- 
based models used in Figures 4-6 and 8-10. Figures 
1-3 were constructed with the same parameters, ex- 
cept for Iq and a. 

model includes two inward Ca 2+ currents (7i and 
Is), an outward delayed-rectifier Potassium current 
(7k), and a leakage current 



C 



.dV 
dT 
dn 
~dt 
ds 
dt 
dc 



= A 



Ic*{V,s) - I K {V,n) - I^{V) 



Tn(V) 
Soo{V, C) - S 



T S (V,C) 

f[-aI Ca (V,s)-k c c] 



(13) 



where V is the membrane potential, c is the intra- 
cellular free Ca 2+ concentration, n and s are the 
activation variables. The ionic currents are given 
by 



h(V)+I s (V,s) 
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h(V) 
h(V,s) 
lK(V,n) 



9i mooOO (V - 
9s s (V - V Ca ) 
g K n (V - Vk) 
9l (Y-Vl) 



V Ca ) 



(14) 



where the sum is taken over all the times to where 
bursts of n spikes are initiated, N n is the total num- 
ber of bursts witn n spikes, and so{t) is the n-burst- 
triggered average (BTA): 



The steady-state functions and activation times are 

1 



s (t) 



to 



to)- 



Xoo(V) 
Soc{V,c) 
Tn(V) 
T S (V,C) 

A{V,c) 



1 



exp [{V x 
1 



V)/S x ] 



for x = m, rtln Equation 16, C p is the N x N 
matrix, 



prior covariancc 



exp[2A(V,c)Y 

1~n0 



1 



exp [{V 

T S 



V n ) /S n ] 



2cosh(A(V,c)) 
V s + S s \og(c)- 



V 



2S S 



and c = c/(l/xM). 



Covariance analysis 



In order to disclose the stimulus features re- 
sponsible for burst generation, we use spike- 
triggered covariance techniques (Paninski, 2003; 
Samengo and Gollisch, 2013). We define P„[io|s(i)] 
as the probability of generating a burst of n spikes 
at time to conditional to a time-dependent stim- 
ulus s(t). We assume that P n only depends on 
the stimulus s(t) through a few relevant features 
f l {t - t Q ), f(t - t ), . . . f k (t - t Q ). The stimulus 
and the relevant features are continuous functions 
of time. For computational purposes, we represent 
them as vectors s and of N components, where 
each component Sj = s(j St) and fj = f % {j St) is 
the value of the stimulus evaluated at discrete inter- 
vals St. If St is small compared with the relevant 
timescales of the models, the discretized signal is 
still a good approximation of the original contin- 
uous function. The relevant features f 1 , . . . , f k lie 
in the space spanned by those eigenvectors of the 
matrix 



(C p )y = s(U + t) s(tj +t)- s(t) 



where the horizontal line represents a temporal av- 
erage on the variable t. 

The eigenvalues of M that are larger than 1 cor- 
^25) respond to directions in stimulus space where the 
stimulus segments associated with burst generation 
have an increased variance, as compared to the raw 
collection of stimulus segments. Such directions are 
often found in cells excited by stimuli that are large 
in modulus, and whose sign is irrelevant. Eigen- 
values that lie significantly below unity are associ- 
ated with stimulus directions of decreased variance. 
Such stimuli are often found in cells that respond 
to precisely tuned values of the stimulus. In general 
terms, the obtained eigenvalue is a measure of the 
ratio of variances of the two ensembles. That is, 
an eigenvalue that is noticeably smaller than unity 
indicates a certain feature for which the ratio of 
the variance of the burst-triggering stimuli and the 
variance of the raw stimuli is significantly small. 
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T 
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Table 3: Parameters of the Chay-Cook model models used in Figures. 4-6 and 8-10. 
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